%% Meager (2019) alternative grouping analysis
% Replicates Table 4

% Random seed
master_seed = 11;
rng(master_seed)
path_num = 10; % number of permutations
seed_list = randi([1,1e6], path_num, 1);

% Initialization
N = 15;
max_iter = 100; 
sim_cov = 5000;
cov_stat_sim = simulate_cov_stat_distribution(N, sim_cov);

EB_nu_list = zeros(path_num, 1);
ecf_benchmark = zeros(path_num, 1);
ecf_p_benchmark = zeros(path_num, 1);
ecf_EB = zeros(path_num, 1);
ecf_p_EB = zeros(path_num, 1);

for k = 1:path_num
    % Alternative grouping
    % Keep the indices that has data
    profit_indices = [1, 2, 3, 4, 5, 6, 7];
    revenue_indices = [1, 2, 3, 4, 5, 6, 7];
    expend_indices = [1, 2, 3, 4, 5, 6, 7];
    consump_indices = [1, 2, 3, 4, 5];
    durables_indices = [2, 3, 4, 5];
    tempt_indices = [1, 2, 3, 4, 5];
    
    % Randomly divide studies into groups
    group_profit = random_pair(profit_indices);
    group_revenue = random_pair(revenue_indices);
    group_expend = random_pair(expend_indices);
    group_consump = random_pair(consump_indices);
    group_durables = random_pair(durables_indices);
    group_tempt = random_pair(tempt_indices);

    % Extract estimates and SE 
    [estimates, se] = extract_est(group_profit, group_revenue, group_expend, ...
                     group_consump, group_durables, group_tempt);
    
    % Randomly assign baseline and validation studies
    [hat_theta, hat_vartheta, rep_base_se2, rep_val_se2] = random_assignment(estimates, se, seed_list(k));
    
    % ECF under nu=0
    [bar_tau_1, bar_V_tau_1] = compute_posterior_moments(hat_theta, rep_base_se2, 0);
    [predict_se_1, empirical_cov_1, cov_p_value_1] = check_CI(hat_vartheta, rep_val_se2, bar_tau_1, bar_V_tau_1, 0, cov_stat_sim);
    ecf_benchmark(k,:) = empirical_cov_1;
    ecf_p_benchmark(k,:) = cov_p_value_1;

    % Estimate nu
    EB_est = MLE_nu(hat_theta, rep_base_se2, 10^-6);
    EB_nu_list(k,:) = EB_est;

    % ECF under estimated nu
    [bar_tau_2, bar_V_tau_2] = compute_posterior_moments(hat_theta, rep_base_se2, EB_est);
    [predict_se_2, empirical_cov_2, cov_p_value_2] = check_CI(hat_vartheta, rep_val_se2, bar_tau_2, bar_V_tau_2, EB_est, cov_stat_sim);
    ecf_EB(k,:) = empirical_cov_2;
    ecf_p_EB(k,:) = cov_p_value_2;


end



